# updated 1/20/2015 to work with R 3.1.2 and lme4 1.1-7
# ======================================================
# Now estimate survival from July 2007 to July 2008
# ======================================================

surv=read.csv("July 2007 census and biomass2 edited.csv",header = T,na.strings=c("NA","na"))
str(surv)
#converting variables in "surv" data to numeric/categorical variables
surv$comp = as.factor(surv$comp)
surv$herb = as.factor(surv$herb)
surv$species = as.factor(surv$species)
surv$bolt08 = as.factor(surv$bolt)
# figure out which rows have the ugly "" factor levels
surv[surv$Comp=="",5:10] = NA
# relevel competition factor
surv$Comp = factor(surv$Comp,levels=c("Low","Medium","High"))
surv$Species = factor(surv$Species,levels=c("BT","TT"))
surv$Herb = factor(surv$Herb,levels=c("Ambient","Reduced"))
surv = surv[complete.cases(surv[,5:10]),]

with(surv,table(Species,survived08))

# Check for complete cases of the variables we need 
# no.leaves07, survived08, bolt08
# alot of missing values for no.leaves07 in bull thistle
with(surv,table(Species,is.na(no.leaves07)))
surv = surv[complete.cases(surv[,c("no.leaves07","survived08","bolt08")]),]
with(surv,table(Species,survived08)) 

# for compatibility, split into tall and bull datasets
surv.tall = surv[surv$Species=="TT",]
surv.bull = surv[surv$Species=="BT",]
surv.tall$plot = paste(surv.tall$block,surv.tall$plot,sep=".")
surv.bull$plot = paste(surv.bull$block,surv.bull$plot,sep=".")


###########################################################################################
## Fit survival models for tall thistle
## convert from logodds scale to parameters that can be used in IPM
###########################################################################################
# now predict veg.bio07 for this dataset. 
# assumes that growth models have been fitted, check for appropriate models
if (!exists("tb.fits")) stop("fit growth models first!")

# now using model averaged predictions
nd = surv.tall
nd$no.leaves08  = nd$no.leaves07
nd$bolt08 = factor(0,levels=levels(grow.tall$bolt08))
results = matrix(NA,nrow=nrow(surv.tall),ncol=length(tb.fits))
for (i in 1:length(tb.fits)){
	# generate predictions from a lmer model using the data from inside the model, mostly. 
	fm = tb.fits[[i]]
	results[,i] = predict(fm,newdata=nd,re.form=~(1|block))
}
# model averaged prediction for veg.bio07 parameter
# make sure we have the right weights vector!!
bic = unlist(sapply(tb.fits,BIC))
dbic = bic-min(bic)
wbic = exp(-dbic/2)/sum(exp(-dbic/2))
surv.tall[,"veg.bio07"] = apply(results,1,function(x)sum(x*wbic))

# fit global model with block and plot effects
st.mm0 = glmer(survived08~veg.bio07 + Comp + Herb + 
		veg.bio07:Comp + veg.bio07:Herb + Comp:Herb + 
		(1|block) + (1|plot), family=binomial, data=surv.tall)
# compare model with variation in herbivory among blocks
st.mm1 = glmer(survived08~veg.bio07 + Comp + Herb +
		veg.bio07:Comp + veg.bio07:Herb + Comp:Herb + 
		(Herb|block) + (1|plot), family=binomial, data=surv.tall)
st.mm2 = glmer(survived08~veg.bio07 + Comp + Herb +
		veg.bio07:Comp + veg.bio07:Herb + Comp:Herb + 
		(1|block) + (Herb|plot), family=binomial, data=surv.tall)

st.mm3 = glmer(survived08~veg.bio07 + Comp + Herb +
		veg.bio07:Comp + veg.bio07:Herb + Comp:Herb + 
		(1|block) + (Comp|plot), family=binomial, data=surv.tall)
st.mm4 = glmer(survived08~veg.bio07 + Comp + Herb + 
		veg.bio07:Comp + veg.bio07:Herb + Comp:Herb + 
		(Comp|block) + (1|plot), family=binomial, data=surv.tall)
# worse, try removing random effects
st.mm5 = glmer(survived08~veg.bio07 + Comp + Herb +
		veg.bio07:Comp + veg.bio07:Herb + Comp:Herb + 
		(1|block), family=binomial, data=surv.tall)
# NO! very bad, remove block?
st.mm6 = glmer(survived08~veg.bio07 + Comp + Herb +
		veg.bio07:Comp + veg.bio07:Herb + Comp:Herb + 
		(1|plot), family=binomial, data=surv.tall)
sapply(list(st.mm0,st.mm1,st.mm2,st.mm3,st.mm4,st.mm5,st.mm6),AIC)
# OK, keep both random effects, but simple

# OK try as BIC model selection to do multi-model inference
ts.models = list("~veg.bio07 + Herb + Comp + veg.bio07:Comp + veg.bio07:Herb + Herb:Comp",
		"~veg.bio07 + Herb + Comp + veg.bio07:Comp + veg.bio07:Herb",
		"~veg.bio07 + Herb + Comp + veg.bio07:Comp + Herb:Comp",
		"~veg.bio07 + Herb + Comp + veg.bio07:Herb + Herb:Comp",
		"~veg.bio07 + Herb + Comp + veg.bio07:Herb",
		"~veg.bio07 + Herb + Comp + Herb:Comp",
		"~veg.bio07 + Herb + Comp + veg.bio07:Comp",
		"~veg.bio07 + Herb + Comp",
		"~veg.bio07 + Herb +  veg.bio07:Herb",
		"~veg.bio07 + Herb",
		"~veg.bio07 + Comp + veg.bio07:Comp",
		"~veg.bio07 + Comp",
		"~Herb + Comp + Herb:Comp",
		"~Herb + Comp",
		"~Herb",
		"~Comp",
		"~veg.bio07",
		"~1")

ts.fits = lapply(ts.models,function(ff)glmer(paste("survived08",ff,"+ (1|block) + (1|plot)"),data=surv.tall,family=binomial))
bic = unlist(sapply(ts.fits,BIC))
dbic = bic-min(bic)
wbic = exp(-dbic/2)/sum(exp(-dbic/2))
data.frame(unlist(ts.models),dbic,wbic)[order(dbic),]
# model 15 has 82% of the weight, and model 10 has 16%

# make a plot
nd = data.frame(veg.bio07=rep(seq(-2.5,0.8,0.1),times=6),
		Comp=factor(rep(c("Low","Medium","High"),each=68),levels=levels(surv.tall$Comp)),
		Herb=factor(rep(rep(c("Ambient","Reduced"),3),each=34)))

results = matrix(NA,nrow=nrow(nd),ncol=length(ts.fits))
for (i in 1:length(ts.fits)){
	# generate predictions from a glmer model using the data from inside the model, mostly. 
	fm = ts.fits[[i]]
	results[,i] = predict(fm,newdata=nd,type="response",re.form=~0)
}
# model averaged prediction for veg.bio07 parameter
nd[,"surv"] = apply(results,1,function(x)sum(x*wbic))

xyplot(surv~veg.bio07|Comp+Herb,data=nd,type="l",
panel=function(...){
		panel.xyplot(...)
		pn = which.packet()
		#ltext(0,0.5,pn[1])
		pick = surv.tall$Comp == levels(surv.tall$Comp)[pn[1]] & 
			 surv.tall$Herb == levels(surv.tall$Herb)[pn[2]]
		ticks = surv.tall[pick,c("veg.bio07","survived08")]
		panel.rug(ticks[ticks$survived08==0,1])
		panel.rug(ticks[ticks$survived08==1,1],y=NULL, regular=FALSE)
	},
ylim=c(0,1))

#------------------
coef.names = names(fixef(ts.fits[[1]]))
coef.results = matrix(NA,nrow=length(ts.fits),ncol=length(coef.names))
coef.se = matrix(NA,nrow=length(ts.fits),ncol=length(coef.names))
for (i in 1:length(ts.fits)){
	coef.results[i,] = sapply(coef.names,function(cf)fixef(ts.fits[[i]])[cf],USE.NAMES=FALSE)
	se = sqrt(diag(vcov(ts.fits[[i]])))
	names(se) = names(fixef(ts.fits[[i]]))
	coef.se[i,] = sapply(coef.names,function(cf)se[cf])
}

coef.avg = apply(coef.results,2,function(x)sum(x*wbic,na.rm=TRUE))
t1 = sweep(coef.results,2,coef.avg)
t2 = sqrt(coef.se^2 + t1^2)
coef.avg.se = apply(t2,2,function(x)sum(wbic*x,na.rm=TRUE))

nd.intercept = data.frame(veg.bio07=rep(0,times=6),
		Comp=factor(rep(c("Low","Medium","High"),each=2),levels=levels(surv.tall$Comp)),
		Herb=factor(rep(c("Ambient","Reduced"),3)))
mm.intercept = model.matrix(as.formula(ts.models[[1]]),data=nd.intercept)

nd.slope = data.frame(veg.bio07=rep(1,times=6),
		Comp=factor(rep(c("Low","Medium","High"),each=2),levels=levels(surv.tall$Comp)),
		Herb=factor(rep(c("Ambient","Reduced"),3)))

mm.slope = model.matrix(as.formula(ts.models[[1]]),data=nd.slope)

mm.slope[,c(1,3:5,9:10)] = 0

int.results = array(NA,dim=c(length(ts.fits),nrow(mm.intercept)))
int.se = matrix(NA,nrow=length(ts.fits),ncol=nrow(mm.intercept))
slope.results = array(NA,dim=c(length(ts.fits),nrow(mm.slope)))
slope.se = matrix(NA,nrow=length(ts.fits),ncol=nrow(mm.intercept))

for (i in 1:length(ts.fits)){
	fef = rep(0,times=length(coef.names))
	names(fef) = coef.names
	vcf = matrix(0,nrow=length(coef.names),ncol=length(coef.names))
	rownames(vcf) = coef.names
	colnames(vcf) = coef.names

	fe = fixef(ts.fits[[i]])
	fef[names(fe)] = fe
	vcf[names(fe),names(fe)] = as.matrix(vcov(ts.fits[[i]]))
	int.results[i,] = mm.intercept %*% fef
	slope.results[i,] = mm.slope %*% fef
	for (j in 1:nrow(mm.intercept)){
		pick = which(mm.intercept[j,]>0)
		int.se[i,j] = sqrt(sum(vcf[pick,pick]))
		pick = which(mm.slope[j,]>0)
		slope.se[i,j] = sqrt(sum(vcf[pick,pick])) # make it the se not the variance
	}
	
}

int.avg = apply(int.results,2,function(x)sum(x*wbic,na.rm=TRUE))
t1 = sweep(int.results,2,int.avg)
t2 = sqrt(int.se^2 + t1^2)
int.avg.se = apply(t2,2,function(x)sum(wbic*x,na.rm=TRUE))

slope.avg = apply(slope.results,2,function(x)sum(x*wbic,na.rm=TRUE))
t1 = sweep(slope.results,2,slope.avg)
t2 = sqrt(slope.se^2 + t1^2)
slope.avg.se = apply(t2,2,function(x)sum(wbic*x,na.rm=TRUE))

for (i in 1:6) {IPM.parameters$surv.intercept[IPM.parameters$species=="TT" & IPM.parameters$comp==nd.intercept[i,"Comp"] & IPM.parameters$herb==nd.intercept[i,"Herb"]]=int.avg[i]}
for (i in 1:6) {IPM.parameters$surv.intercept.se[IPM.parameters$species=="TT" & IPM.parameters$comp==nd.intercept[i,"Comp"] & IPM.parameters$herb==nd.intercept[i,"Herb"]]=int.avg.se[i]}
for (i in 1:6) {IPM.parameters$surv.slope[IPM.parameters$species=="TT" & IPM.parameters$comp==nd.slope[i,"Comp"] & IPM.parameters$herb==nd.intercept[i,"Herb"]]=slope.avg[i]}
for (i in 1:6) {IPM.parameters$surv.slope.se[IPM.parameters$species=="TT" & IPM.parameters$comp==nd.slope[i,"Comp"] & IPM.parameters$herb==nd.intercept[i,"Herb"]]=slope.avg.se[i]}



###########################################################################################
## Fit survival models for Bull thistle
###########################################################################################
# now predict veg.bio07 for this dataset. 
# assumes that growth models have been fitted, check for appropriate models
if (!exists("bb.fits")) stop("fit growth models first!")

# now using model averaged predictions
nd = surv.bull
nd$no.leaves08  = nd$no.leaves07
nd$bolt08 = factor(0,levels=levels(grow.bull$bolt08))
results = matrix(NA,nrow=nrow(surv.bull),ncol=length(bb.fits))
for (i in 1:length(bb.fits)){
	# generate predictions from a lmer model using the data from inside the model, mostly. 
	fm = bb.fits[[i]]
	results[,i] = predict(fm,newdata=nd,re.form=~(1|block))
}
# model averaged prediction for veg.bio07 parameter
# make sure we have the right weights
bic = unlist(sapply(bb.fits,BIC))
dbic = bic-min(bic)
wbic = exp(-dbic/2)/sum(exp(-dbic/2))
surv.bull[,"veg.bio07"] = apply(results,1,function(x)sum(x*wbic))

# fit global model with block and plot effects
sb.mm0 = glmer(survived08~veg.bio07 + Comp + Herb + 
		veg.bio07:Comp + veg.bio07:Herb + Comp:Herb + 
		(1|block) + (1|plot), family=binomial, data=surv.bull)
# compare model with variation in herbivory among blocks
sb.mm1 = glmer(survived08~veg.bio07 + Comp + Herb +
		veg.bio07:Comp + veg.bio07:Herb + Comp:Herb + 
		(Herb|block) + (1|plot), family=binomial, data=surv.bull)
sb.mm2 = glmer(survived08~veg.bio07 + Comp + Herb +
		veg.bio07:Comp + veg.bio07:Herb + Comp:Herb + 
		(1|block) + (Herb|plot), family=binomial, data=surv.bull)

sb.mm3 = glmer(survived08~veg.bio07 + Comp + Herb +
		veg.bio07:Comp + veg.bio07:Herb + Comp:Herb + 
		(1|block) + (Comp|plot), family=binomial, data=surv.bull)
sb.mm4 = glmer(survived08~veg.bio07 + Comp + Herb + 
		veg.bio07:Comp + veg.bio07:Herb + Comp:Herb + 
		(Comp|block) + (1|plot), family=binomial, data=surv.bull)
# worse, try removing random effects
sb.mm5 = glmer(survived08~veg.bio07 + Comp + Herb +
		veg.bio07:Comp + veg.bio07:Herb + Comp:Herb + 
		(1|block), family=binomial, data=surv.bull)
# NO! very bad, remove block?
sb.mm6 = glmer(survived08~veg.bio07 + Comp + Herb +
		veg.bio07:Comp + veg.bio07:Herb + Comp:Herb + 
		(1|plot), family=binomial, data=surv.bull)
sapply(list(sb.mm0,sb.mm1,sb.mm2,sb.mm3,sb.mm4,sb.mm5,sb.mm6),AIC)
# OK, keep both random effects, but simple

# OK try as BIC model selection to do multi-model inference
bs.models = list("~veg.bio07 + Herb + Comp + veg.bio07:Comp + veg.bio07:Herb + Herb:Comp",
		"~veg.bio07 + Herb + Comp + veg.bio07:Comp + veg.bio07:Herb",
		"~veg.bio07 + Herb + Comp + veg.bio07:Comp + Herb:Comp",
		"~veg.bio07 + Herb + Comp + veg.bio07:Herb + Herb:Comp",
		"~veg.bio07 + Herb + Comp + veg.bio07:Herb",
		"~veg.bio07 + Herb + Comp + Herb:Comp",
		"~veg.bio07 + Herb + Comp + veg.bio07:Comp",
		"~veg.bio07 + Herb + Comp",
		"~veg.bio07 + Herb +  veg.bio07:Herb",
		"~veg.bio07 + Herb",
		"~veg.bio07 + Comp + veg.bio07:Comp",
		"~veg.bio07 + Comp",
		"~Herb + Comp + Herb:Comp",
		"~Herb + Comp",
		"~Herb",
		"~Comp",
		"~veg.bio07",
		"~1")

bs.fits = lapply(bs.models,function(ff)glmer(paste("survived08",ff,"+ (1|block) + (1|plot)"),data=surv.bull,family=binomial))
bic = unlist(sapply(bs.fits,BIC))
dbic = bic-min(bic)
wbic = exp(-dbic/2)/sum(exp(-dbic/2))
data.frame(unlist(bs.models),dbic,wbic)[order(dbic),]
# model 10 has 83% of the weight, and model 8 has 11%

# make a plot
nd = data.frame(veg.bio07=rep(seq(-2.5,0.8,0.1),times=6),
		Comp=factor(rep(c("Low","Medium","High"),each=68),levels=levels(surv.bull$Comp)),
		Herb=factor(rep(rep(c("Ambient","Reduced"),3),each=34)))

results = matrix(NA,nrow=nrow(nd),ncol=length(bs.fits))
for (i in 1:length(bs.fits)){
	# generate predictions from a glmer model using the data from inside the model, mostly. 
	fm = bs.fits[[i]]
	results[,i] = predict(fm,newdata=nd,re.form=~0,type="response")
}
# model averaged prediction for veg.bio07 parameter
nd[,"surv"] = apply(results,1,function(x)sum(x*wbic))

xyplot(surv~veg.bio07|Comp+Herb,data=nd,type="l",
panel=function(...){
		panel.xyplot(...)
		pn = which.packet()
		#ltext(0,0.5,pn[1])
		pick = surv.bull$Comp == levels(surv.bull$Comp)[pn[1]] & 
			 surv.bull$Herb == levels(surv.bull$Herb)[pn[2]]
		ticks = surv.bull[pick,c("veg.bio07","survived08")]
		panel.rug(ticks[ticks$survived08==0,1])
		panel.rug(ticks[ticks$survived08==1,1],y=NULL, regular=FALSE)
	},
ylim=c(0,1))


#------------------
coef.names = names(fixef(bs.fits[[1]]))
coef.results = matrix(NA,nrow=length(bs.fits),ncol=length(coef.names))
coef.se = matrix(NA,nrow=length(bs.fits),ncol=length(coef.names))
for (i in 1:length(bs.fits)){
	coef.results[i,] = sapply(coef.names,function(cf)fixef(bs.fits[[i]])[cf],USE.NAMES=FALSE)
	se = sqrt(diag(vcov(bs.fits[[i]])))
	names(se) = names(fixef(bs.fits[[i]]))
	coef.se[i,] = sapply(coef.names,function(cf)se[cf])
}

coef.avg = apply(coef.results,2,function(x)sum(x*wbic,na.rm=TRUE))
t1 = sweep(coef.results,2,coef.avg)
t2 = sqrt(coef.se^2 + t1^2)
coef.avg.se = apply(t2,2,function(x)sum(wbic*x,na.rm=TRUE))

nd.intercept = data.frame(veg.bio07=rep(0,times=6),
		Comp=factor(rep(c("Low","Medium","High"),each=2),levels=levels(surv.tall$Comp)),
		Herb=factor(rep(c("Ambient","Reduced"),3)))
mm.intercept = model.matrix(as.formula(bs.models[[1]]),data=nd.intercept)

nd.slope = data.frame(veg.bio07=rep(1,times=6),
		Comp=factor(rep(c("Low","Medium","High"),each=2),levels=levels(surv.tall$Comp)),
		Herb=factor(rep(c("Ambient","Reduced"),3)))

mm.slope = model.matrix(as.formula(bs.models[[1]]),data=nd.slope)

mm.slope[,c(1,3:5,9:10)] = 0



int.results = array(NA,dim=c(length(bs.fits),nrow(mm.intercept)))
int.se = matrix(NA,nrow=length(bs.fits),ncol=nrow(mm.intercept))
slope.results = array(NA,dim=c(length(bs.fits),nrow(mm.slope)))
slope.se = matrix(NA,nrow=length(bs.fits),ncol=nrow(mm.intercept))

for (i in 1:length(bs.fits)){
	fef = rep(0,times=length(coef.names))
	names(fef) = coef.names
	vcf = matrix(0,nrow=length(coef.names),ncol=length(coef.names))
	rownames(vcf) = coef.names
	colnames(vcf) = coef.names

	fe = fixef(bs.fits[[i]])
	fef[names(fe)] = fe
	vcf[names(fe),names(fe)] = as.matrix(vcov(bs.fits[[i]]))
	int.results[i,] = mm.intercept %*% fef
	slope.results[i,] = mm.slope %*% fef
	for (j in 1:nrow(mm.intercept)){
		pick = which(mm.intercept[j,]>0)
		int.se[i,j] = sqrt(sum(vcf[pick,pick]))
		pick = which(mm.slope[j,]>0)
		slope.se[i,j] = sqrt(sum(vcf[pick,pick])) # make it the se not the variance
	}
	
}

int.avg = apply(int.results,2,function(x)sum(x*wbic,na.rm=TRUE))
t1 = sweep(int.results,2,int.avg)
t2 = sqrt(int.se^2 + t1^2)
int.avg.se = apply(t2,2,function(x)sum(wbic*x,na.rm=TRUE))

slope.avg = apply(slope.results,2,function(x)sum(x*wbic,na.rm=TRUE))
t1 = sweep(slope.results,2,slope.avg)
t2 = sqrt(slope.se^2 + t1^2)
slope.avg.se = apply(t2,2,function(x)sum(wbic*x,na.rm=TRUE))



for (i in 1:6) {IPM.parameters$surv.intercept[IPM.parameters$species=="BT" & IPM.parameters$comp==nd.intercept[i,"Comp"] & IPM.parameters$herb==nd.intercept[i,"Herb"]]=int.avg[i]}
for (i in 1:6) {IPM.parameters$surv.intercept.se[IPM.parameters$species=="BT" & IPM.parameters$comp==nd.intercept[i,"Comp"] & IPM.parameters$herb==nd.intercept[i,"Herb"]]=int.avg.se[i]}
for (i in 1:6) {IPM.parameters$surv.slope[IPM.parameters$species=="BT" & IPM.parameters$comp==nd.slope[i,"Comp"] & IPM.parameters$herb==nd.intercept[i,"Herb"]]=slope.avg[i]}
for (i in 1:6) {IPM.parameters$surv.slope.se[IPM.parameters$species=="BT" & IPM.parameters$comp==nd.slope[i,"Comp"] & IPM.parameters$herb==nd.intercept[i,"Herb"]]=slope.avg.se[i]}

write.csv(IPM.parameters,"IPM.Parameters.csv",row.names = FALSE)